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Abstract 

Interior  point  algorithms  for  solving  linear  programming  problems  are  considered.  An 
optimal  3-dimensional  method  is  developed  that,  at  each  iteration,  solves  a subproblem 
bcised  on  minimizing  the  cost  function  on  low-dimensional  cross  sections  of  the  feasible 
region.  This  idea  was  used  by  the  authors  in  [Boggs  et  al.,  Algorithmic  Enhancements  to 
the  Method  of  Centers  for  Linear  Programming  Problems,  ORSA  Journal  of  Computing, 
l(3):159-171,  1989.]  The  generators  for  the  original  2-dimensional  subproblems  are  derived 
from  either  a discrete  or  a continuous  version  of  Huard’s  method  of  centers.  The  generators 
for  the  optimal  3-dimensional  subproblem  include  the  dual  affine  search  direction,  and  two 
higher-order  search  directions.  One  of  the  higher  order  directions  is  a third  order  correction 
to  the  Newton  recentering  direction,  and  the  other  is  a correction  to  the  dual  affine  direction 
that  is  motivated  by  the  use  of  rank-one  updates  of  the  second  derivative  information. 
Numerical  results  are  presented  for  the  optimal  3-dimensional  method  that  indicate  almost 
a 18.9%  reduction  in  CPU  time  compared  to  our  best  dual  affine  implementation. 
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1.  Introduction 


In  our  previous  work  on  solving  linear  programming  (LP)  problems  [BogDDW89],  we  suggested 
an  interior  point  strategy  that  consists  of  exactly  solving  a sequence  of  2-dimensional  subprob- 
lems. These  2-dimensional  subproblems  are  generated  by  considering  the  restriction  of  the 
original  LP  to  two  search  directions.  In  our  best  2-dimensional  algorithm,  one  of  the  directions 
is  the  dual  affine  direction  and  the  second  direction  is  an  approximate  recentering  direction.  In 
this  paper,  we  extend  this  idea  to  include  a third,  higher-order  direction  derived  from  Huard’s 
method  of  centers  [Hua67].  Our  numerical  results  demonstrate  the  effectiveness  of  using  such 
higher-order  information:  the  number  of  iterations  required  to  solve  a subset  of  the  Netlib  prob- 
lems [Gay85]  (those  problems  without  explicit  bounds)  is  reduced  by  33.4%  and  the  time  by 
18.9%  over  our  new  version  of  the  dual  affine  method.  Our  results  are  competitive  with  the 
dual  affine  procedures  reported  in  [Meh89],  [MonM87]  and  [AdlRV86],  and  with  the  primal-dual 
interior  point  methods  reported  in  [Sha85],  [McSMS88]  and  [LusMS89]. 

Our  techniques,  which  we  call  optimal  multi-dimensional  methods,  effectively  solve  the  prob- 
lem of  how  to  combine  the  various  search  directions  that  arise  from  interior  point  methods.  The 
theoretical  and  algorithmic  details  of  our  optimal  3-dimensional  methods  are  contained  in  §2, 
and  the  implementation  details  are  given  in  §3.  The  computational  results  for  the  subset  of  the 
Netlib  test  problems  [Gay85]  are  presented  in  §4. 

In  the  remainder  of  this  section,  we  motivate  our  optimal  2-dimensional  method  using 
Huard’s  original  method  of  centers  and  a continuous  version  of  it.  The  various  solution  “tra- 
jectories” discussed  below  are  covered  in  greater  detail  in  [BayL86]  and  [WitBD88b].  The 
corresponding  numerical  procedures  evolving  from  them  are  discussed  in  [BogDDW89]. 

The  LP  problem  that  we  solve  is  of  the  form 


T 

mm  c u 

XL 

s.t.  Au  < b 


(1.1) 


where  c,  u G A G and  h G R’”.  We  define  the  set  of  residuals  related  to  the  con- 

straints of  ( 1.1)  as 

rfc(u)  = bk  - AkU,  k - 1, . . . , m, 

where  Ak  denotes  the  kth  row  of  A,  and  let  ro(u,  t)  = t — correspond  to  the  residual  of  the 
objective  row  assuming  an  upper  bound  t on  the  objective  function  value.  In  particular,  let  uq 
be  a feasible  point.  Then  for  to  = c'^uq  + e,  e > 0, 


ro{uo,to)>0  and  (uo)  > 0,  A:  = 1, . . . , m. 
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The  center,  defined  by  the  constraints  of  (1.1)  and  the  objective  constraint  for  t = to,  is  the 
feasible  point  u that  solves 

i^(to)  = max[L(u,  to)l  (1-2) 


where 


L{u,to)  = logro(u,to)  + ^ log  rfc(u). 


(1.3) 


The  optimization  problem  (1.2)  can  be  solved  using  Newton’s  method.  The  search  direction  is 
then  given  by 

Sn  = [Vu„i(u,  t)]“W„L(u,  t) 

T1-1 (1-4) 


A'^D^A  + ^ 


1 1 


, R = (^,  . . . , and  e = t — c”u  > 0.  Applying  the 


where  D = diag 
Sherman-Morrison  formula  to  (1.4)  results  in 

= /?i(e)  [A^ A]~\  + ^2{e)  [a'^ D'^ AY^  A^ R 


:i.5) 


for  real-valued  functions  /3i(e)  and  y92(€).  The  first  term  in  (1.5)  is  the  dual  affine  search  direction 
[AdlRV86]  and  the  second  is  the  Newton  search  direction  for  locating  the  center  of  the  original 
polytope  defined  by  Au  < h without  the  objective  constraint. 

Note  that  one  could  systematically  reduce  the  value  of  t and  solve  a sequence  of  centering 
problems  (1.2).  This  yields  a sequence  of  iterates  Ui  with  successively  lower  objective  function 
values.  With  a reasonable  selection  of  t,,  it  can  be  shown  that  {tq}  converges  to  an  optimal 
solution  as  z — ♦ oo.  This  procedure  is  Huard’s  original  method  of  centers  [Hua67]  applied  to 
the  linear  programming  problem.  The  implementation  of  Huard’s  method  proposed  by  Renegar 
[Ren86]  was  shown  to  possess  an  equivalent  polynomial  complexity  bound  to  that  of  Karmarkar’s 
original  method  [Kar84]. 

It  is  also  easily  seen  that  by  continuously  moving  the  constraint  corresponding  to  the  objective 
function,  one  can  find  a continuous  trajectory,  rather  than  a discrete  set  of  points  Ufc,  that 
converges  to  an  optimal  solution  of  (l.l).  In  particular,  differentiating  (1.4)  with  respect  to  t 
results  in  the  ODE 

u[t]  ~ -VuuL(u,  t)"^VutI(u,  t).  (1.6) 

This  ODE  can  be  supplied  with  initial  conditions  consisting  of  an  arbitrary  feasible  point  and 
an  appropriate  value  of  t with  the  result  that  for  every  feasible  point  there  is  a trajectory  that 
connects  it  to  an  optimal  solution  of  (l.l).  The  direction  field  at  each  feasible  point  can  be 
easily  shown  to  be  proportional  to  the  dual  affine  direction  described  earlier. 
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Since  the  solutions  of  (1.6)  can  get  arbitr<irily  close  to  an  exponential  number  of  vertices 
(see,  e.g.,  [MegS86]),  it  is  of  interest  to  consider  modifications  that  “correct”  the  trajectories 
towards  the  center  of  the  polytope.  Thus  a second  ODE  can  be  formulated  [BogDDW89]  that 
includes  the  Newton  recentering  component  (1.5),  i.e., 

u'(t)  = - Vuul,(u,  t)~^  [Vutl,(u,  t)  - 0VuL(u,  t)]  (1.7) 

for  (^  > 0.  If  VuL{u,t)  = 0,  then  the  recentering  component  has  no  effect,  and  the  resulting 
solution  is  referred  to  as  the  center  trajectory. 

Both  (1.5)  and  (1.7)  indicate  the  motivation  for  combining  the  two  directions 

31  = {A^D^Ay\  (1.8) 

32  = A^R  (1.9) 

into  a single  search  direction 

3 = ClSl  + ^252- 

The  values  of  the  weights  and  f2  dramatically  affect  the  performance  of  an  interior  point 
procedure.  One  successful  method  for  determining  how  to  combine  such  multiple  search  direc- 
tions when  using  large  steplengths  is  the  subproblem  approach.  This  approach  was  suggested 
by  [Kar85]  and  [Gon87a],  and  explored  more  fully  in  [BogDDW89]. 

Initially,  we  considered  the  2-dimensional  subproblem  defined  by  the  span  of  the  two  search 
directions  (1.8)  and  (1.9).  Assuming  these  directions  are  not  colinear,  i.e.,  the  current  point  u is 
feasible  and  not  on  the  center  trajectory,  3i  and  S2  define  a 2-dimensional  plane  that  intersects 
the  fuU  LP  polytope.  The  full  LP  objective  c^u  can  then  be  minimized  on  this  plane  using  the 
linear  program; 

min  ^ic^3i  + ^23^ 32 
subject  to 

^lAsi  + f2>f52  <h  — An 
> 0 

^2  unrestricted. 

Note  that  since  u is  assumed  to  be  feasible,  = 0 and  ^2  = 0 is  also  feasible.  Further,  since 
A^D^A  is  positive  definite,  the  dual  affine  direction  3i  will  always  be  a descent  direction  for  the 
objective  function;  hence,  the  sign  restriction  on  These  two  conditions  imply  that  a cost- 
improving feasible  solution  always  exists  for  the  subproblem.  Finally,  observe  that  the  resulting 
weights  and  ^2  are  not  dependent  on  the  arbitrary  values  of  e and  <f>  defined  earlier. 
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2.  New  Algorithms 

2.1.  Motivation 

The  optimal  multi-dimensional  approach  is  motivated  by  noting  that  the  cost  function  is  eas- 
ily mmimked  on  a low-dimensional  polytope,  and  that  the  solution  to  this  reduced  problem 
determines  a search  direction  that  optimally  weights  the  component  search  directions  over  the 
low-dimensional  cross  section  of  the  full  LP  poly  tope.  The  new  algorithm  presented  in  this 
section  uses  a 3-dimensional  subproblem.  The  new  search  direction  included  in  this  subproblem 
can  be  derived  either  from  a third-order  correction  to  the  Newton  recentering  direction  (1.5),  or 
from  a higher-order  technique  for  solving  (1.7). 

2.2.  The  Subproblem  Generators 

The  sole  restriction  on  the  subproblem  generators  is  that  they  be  linearly  independent.  In  (1.10), 
the  two  generators 

31  = {A^D^Ay\  (2.1) 

32  = {A^D^Ay^  A'^R.  (2.2) 

are  linearly  independent  unless  the  current  estimate  u lies  on  the  center  trajectory,  in  which 
case,  the  dual  affine  direction  (2.1)  is  the  obvious  search  direction.  In  [BogDDW89],  we  discuss 
other  possibilities  for  32,  including  the  first-order  (steepest  descent)  recentering  direction 

32  = A'^R. 


We  also  derive  a correction  to  the  dual  affine  direction, 

S2=  {A’^D^Ay^  aJ,  (2.3) 

where  k is  the  index  of  the  first  constraint  encountered  in  the  dual  affine  direction.  This  direction 
is  motivated  by  considering  the  search  direction  3i  = (A^D'^A)  ^ c,  and  the  rank-one  update 
of  [A'^D^A)  ^ reflecting  the  change  in  residual  k which  produces  the  “new”  search  direction 
^new  _ D‘2 j^y^^  g.  Assuming  c is  not  orthogonal  to  Ak,  it  is  easily  shown  that 

3""“'  = 0131  + 02  {A'^D^Ay'-  AJ 

for  scalars  Oi  and  02.  Thus,  the  direction  [A'^ D^A)  ^ A'^  represents  a rank-one  correction  to 
the  dual  affine  direction. 
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There  cire  several  ways  to  derive  the  third  direction,  S3,  for  our  optimal  3-dimensional  method. 
By  recalling  the  differential  equation  (1.7),  one  can  expand  the  solution  u{t)  in  a Taylor  series. 
The  directions  for  the  optimal  2-dimensional  method  can  then  be  viewed  as  the  independent 
directions  that  arise  from  considering  the  expansion  through  two  terms.  By  considering  the 
expansion  through  three  terms,  the  new  direction  S3  = s^,  described  in  detail  in  §2.3  below, 
is  obtained.  It  is  easily  shown  that  this  new  direction  can  be  found  by  considering  either  the 
expansion  based  on  the  dual  affine  direction  or  the  Newton  recentering  component.  Since  the 
latter  is  easier  to  develop,  we  present  that  derivation  in  §2.3.  We  also  show  in  §2.4  that  this 
direction  can  be  economically  calculated  by  using  the  “factorable”  form  of  the  terms  involved. 

The  set  of  generators  for  the  subproblem  can  be  varied  to  create  a globally  effective  algorithm. 
See  §3  for  the  details  of  our  implementation  of  this  algorithm,  and  [BogDDW89j  for  a detailed 
description  and  earlier  computational  results  using  the  2-dimensional  subproblem  generators. 


2.3.  Third  Order  Correction  Terms 


The  derivation  of  the  term  s/i  as  a correction  to  the  recentering  direction  follows  that  described 
in  [JacM86].  To  simplify  the  presentation,  we  give  an  informal  derivation  in  one  dimension;  the 
extension  to  higher  dimensions  is  straightforward.  Thus,  assume  that  we  fix  t and  that  we  use 
a Taylor’s  series  expansion  of  the  scalar  function  L[u)  = L(u,  t), 

L[u  -h  8u)  = L[u)  -f  L'{u)6u  + L"[u)^—^ — t-  L"'[u)  ^ ^ — h . . . . (2.4) 


The  standard  second-order  Newton  method  for  maximizing  the  function  L[u)  assumes  that 

((5u)3 


V"[u)- 


3! 


0, 


resulting  in 

»„  = |L"Hr‘L'(u). 

This  direction  Sn  thus  maximizes  (2.4)  through  three  terms.  Instead,  suppose  that  a significant 
third-order  correction  term  exists,  say  Sh-  Then  for  6u  = Sn  + Sh,  the  recentering  optimality 
condition  requires 

L'[u  + 6u)  = 0, 


33  before.  Ignoring  terms  of  order  Odl^uH"*)  in  (2.4),  we  have 

0 = L'M  + L''[u){3„  + 3U)  + £"'(u)  + 0(||faf ) 

= L'M  + r (u)(.„  + 3,)  + + 0(||«uf ) 
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Note  that  s„  is  C)(|jL'(u)||).  Assuming  sh  is  0(||L'(«)||^),  then  SnSh  is  0(||i/'(u)||^)  and  is 
0(||X'(u)||^).  Ignoring  C>(|lL'(u)||^)  terms,  we  have 

0 = I'(u)  + L"(u)(s„  + s^)  + r»(«)^ 

= V{u)-L'{u)  + L"{u)3^  + L^'’{u)^. 

Hence,  the  third-order  correction  term  is 

= -i|L"(„)|-'L"'(u),^  (2.5) 

2.4.  The  Practorable  Form  of  L''*{u)3^ 

The  actual  computation  of  3h  is  considerably  simplified  using  the  outer  product  form  of  L'". 
With  the  notation  from  [McCS88],  we  consider 


g{x)  = U{a{x)), 


for  U[a)  possessing  as  many  continuous  derivatives  as  required  and  a = a^x.  Then 


Vg(x)^3  = [8U[a[x])/6a.\a^  3 

V^gr(i)s  = a\8^U[oi[x))  / 8a^\[a^  s) 

V^g(i)  (8  3^  = a[<i^l7(a{i))/(5Q:^](a^3)^. 

Let  U{ak{u))  = log(6fc  — Aku).  Then  the  third-order  correction  term  (2.5)  can  be  written  as 


1 


k=l 


{bk  - Afcu)3 


{AkSr 


(2.6) 


The  third-order  correction  term,  3^,  thus  requires  no  increase  in  the  order  of  work  per  iteration 
since  the  Cholesky  factorization  of  L"[u)  has  been  previously  computed  and  since  the  summation 
in  (2.6)  is  0(nm)  operations. 


2.5.  Solving  the  Subproblem. 

The  subproblem  solution  is  found  using  a specialized  revised  dual  simplex  approach  applied 
to  the  ducil  formulation  of  the  subproblem.  This  procedure  is  used  to  minimize  the  dimension 
of  the  basis  matrix  and  to  provide  feasible  primal  solutions  to  the  original  subproblem.  The 
procedure  also  aUows  the  user  to  terminate  the  subproblem  solution  process  after  a fixed  number 
of  iterations  with  a feasible  set  of  weights  for  the  search  direction. 
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The  actual  primal  subproblem  consists  of  five  columns  cind  m rows,  (m  » 5).  The  five 
columns  correspond  to  the  positive  dual  aifine  direction  (2.1),  the  third-order  correction  direction 
(2.5),  cind  either  the  recentering  direction  (2.2)  or  the  rank-one  update  direction  (2.3)  described 
earlier.  The  second  and  third  columns  are  replicated  and  negated  in  the  subproblem  to  allow 
for  positive  and  negative  multipliers.  Since  the  zero  vector  is  feasible  to  the  subproblem,  the 
primal  problem  is  always  initially  feasible. 

The  revised  dual  simplex  approach  is  applied  to  the  dual  formulation  of  the  subproblem. 
Let  S denote  the  n x 5 matrix  composed  of  search  directions  Si, . . . , 35.  Then  dual  subproblem 
formulation  is 

T 

mm  r y 
y 

subject  to  ^2  'j\ 

— [AS)'^y  < S'^ c 
y >0, 

where  the  zth  constraint  of  (2.7)  is  (— < sjc,  for  i = 1,  ...,5.  This  dual  LP  is  dual 
feasible  since  ry  > 0,  for  all  j. 

The  procedure  begins  by  forming  the  LU  factorization  of  the  5x5  basis  matrix  B (initially 
B = I)  and  finding  the  updated  solution  vector  B~^S'^c.  Next,  the  dual  variables,  0^  = rBB-\ 
are  computed  where  is  the  row  vector  of  objective  coefficients  corresponding  to  columns  in 
the  current  basis.  The  pivot  row  selection  rule  finds  the  row  with  the  largest  index  k such  that 
5-5^3  < 0,  corresponding  to  an  infeasible  primal  variable,  where  the  zth  row  of  B ^ is  denoted 
as  B~^.  Because  m » 5,  the  A:th  row  of  B~^  is  computed  explicitly  and  is  used  to  update  the 
m entries  in  row  k.  Computing  B^^  requires  only  three  backsolves  using  the  LU  factorization 
of  B since  two  columns  of  are  known  to  be  unit  vectors. 

The  updated  Hh  row  of  ( — ^45)^,  B^^(  — A^)^,  is  computed  and  column  j is  selected  that 
minimizes 

^T(-^g)T-r, 

over  z such  that  B'^^[  — AS)  I < 0.  where  ( — A5)^  is  the  zth  column  of  ( — A5)^.  This  column 
selection  rule  guarantees  dual  feasibility  of  /?  and  hence  provides  feasible  search  direction  mul- 
tipliers at  each  step.  The  pivot  is  performed  and  the  basis  is  updated  to  include  column  j. 
This  process  is  repeated  until  B~^S'^c  > 0.  For  complete  details  of  the  dual  simplex  approach 
described  above  see  Bazaraa  and  Jarvis  [BcizJ77]. 

Empirical  evidence  suggests  the  subproblem  polytopes  often  have  a large  number  of  redun- 
dcint  constraints.  Using  the  specialized  procedure  described  above,  the  optimal  vertex  is  usually 
found  after  a very  small  number  of  simplex  pivots.  (For  the  test  problems  reported  in  §4,  all 
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subproblem  solutions  required  fewer  than  50  pivots,  and  more  than  95%  of  the  total  number  of 
subproblem  solutions  required  fewer  than  10  pivots.)  An  inordinate  number  of  pivots  could  oc- 
ctir,  however,  due  to  the  redundant  constraints  at  a near-optimal  vertex.  Suboptimal  solutions, 
e.g.,  those  obtained  after  a fixed  number  of  simplex  pivots,  may  be  substituted  for  the  optimal 
solution  to  save  computation  time.  For  the  test  problems  reported  in  §4,  however,  the  sub- 
problem solution  time  is  only  10%  of  the  total  optimal  3-dimensional  solution  time;  terminating 
the  subproblem  prematurely  causes  an  increase  in  the  number  of  major  iterations  resulting  in  a 
overall  increase  in  CPU  time. 

3.  Implementational  Details. 

Optimal  3-Dimensional  Method.  In  the  results  reported  here,  we  use  the  following  basic 
algorithm. 

Algorithm  3.1.  Optimal  3-Dimensional  Method 

1;  Compute  si  = (A^D^A)  at  u^. 

2:  Compute  either 

2.1:  the  rank-one  update  direction  by 

a:  finding  the  index  k of  the  first  constraint  encountered  in  the  direction  si  and 

b:  computing  S2  = (A'^D'^A)  ^ At,  or 

2.2:  the  recentering  direction  S2  = [A'^ D^A)  ^ A^ R. 

3:  Compute  the  third-order  correction  direction  33  = s/i  as  in  (2.6). 

4:  Solve  for  ^1,  ^2  S’Hd  ^3  in  the  3-dimensional  subproblem  defined  by  3i,  32  and  33.  i.e., 

min  {^ic’^3i  + f2c’^S2  + f3c’’^53} 

fl.f3.f3 

subject  to 

$■1  A3i  + ^2^32  + $“3^33  < b — Au'  (^-1) 

>0 

^2,^3  unrestricted. 

5:  Compute  u,-+i  = u,  + (0.99){fi3i  + ^-252  + ^■353]- 

Problem  Scaling.  Our  scaling  algorithm  for  A uses  techniques  discussed  in  detail  in 
[GilMWSl,  p.353].  First,  each  row  of  A is  scaled  by  the  geometric  mean  of  the  minimum 
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and  maximum  absolute  values  of  the  nonzero  elements  of  the  row.  The  columns  of  A are  next 
scaled  in  the  analogous  manner.  This  scaling  of  the  rows  and  columns  is  then  repeated  until  the 
greatest  ratio  of  two  nonzero  elements  in  the  same  column  does  not  change  by  more  than  10%. 
Finally,  each  row  is  scaled  by  its  meiximum  magnitude,  followed  by  a scaling  of  each  column 
by  its  maximum  magnitude.  Although  this  scaling  of  A requires  multiple  passes  through  the 
data  structures,  the  actual  cost  is  minimal  as  is  seen  in  §4.  We  find  this  scaling  improves  the 
robustness  both  of  our  optimal  3-dimensional  method  and  of  our  dual  affine  procedure  over  that 
observed  using  scaling  with  a fixed  number  of  passes  as  recently  suggested  in  [MarSSPBSS]. 

The  subproblem  constraint  matrix  defined  in  (3.1)  is  also  scaled.  First,  each  column  of  the 
subproblem  constraint  matrix  is  scaled  by  the  square  of  the  two-norm  of  the  corresponding 
subproblem  generator  Sk,  k = 1,2,3,  (see  §3).  Then  each  row  is  scaled  by  its  two-norm.  This 
scaling  of  the  subproblem  constraint  matrix  was  also  found  to  improve  the  numerical  stability 
of  the  method. 

Starting  Values  and  Phase  1 Procedure.  Initially,  uq  = Q.  We  observe  that,  for  the 
problems  examined,  this  choice  of  starting  value  outperforms  starting  values  derived  as  a function 
of  the  problem  data.  (See,  e.g.,  [AdlRV86],  [MarSSPB88|,  and  [Meh89].)  Choi  et  al.  [ChoMS88] 
reports  similar  experience  with  uq  = 0.  Lustig  et  al.  [LusMS89],  however,  notes  that  this  choice 
may  not  be  satisfactory  for  the  full  expanded  Netlib  test  set. 

When  necessary,  an  initial  feasible  solution  is  obtained  using  a big-M  Phase  1 procedure  (see, 
e.g.,  [BazJ77]).  We  implement  this  by  appending  a dense  column  to  the  constraint  matrix  A to 
produce  the  Phase  1 problem 

min  c'^u  Mua 

U,Ua 

S.t.  Au  — Uat  < h 

for  e = (1,...,!)^.  (The  initial  choice  for  is  discussed  below.)  As  was  true  in  our  earlier 
work  [BogDDW89],  the  value  of  M is  not  dynamically  updated  at  each  iteration.  Unlike  our 
earlier  work,  however,  the  coefficient  M of  the  artificial  variable  is  a function  of  the  scaled  values 
of  c,  rather  than  always  the  value  10*. 

It  has  been  widely  reported  that  interior  point  methods  are  sensitive  to  the  selection  of  M. 
Numerical  difficulities  cire  encounted  when  M is  either  too  Icirge  or  too  small.  For  the  primal-dual 
procedures,  Choi  et  al.  [ChoMS88]  uses  M = pn^ max{||c||jjjj , ||6||jjq},  where  p is  either  3 or  30. 
Similarly,  Mehrotra  [Meh89]  uses  M = 10^  ||c||.  For  our  work,  we  set 

M = 10  X min{l0^,  ||c||qq  x max{l0*,  ||c||qq}}. 
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We  found  this  heuristic  more  robust  than  those  previously  reported.  It  is  arbitrary,  but  when 
coupled  with  our  scaling  of  A,  its  use  significantly  improves  both  our  optimal  3-dimensional 
results,  and  to  a lesser  degree,  our  dual  affine  results  (see  §4). 

In  addition  to  the  selection  of  M,  the  big-M  Phase  1 procedure  requires  that  we  initialize  the 
artificial  Vciriable  Ua-  In  our  earlier  work,  we  set  Ua  = 2 ||r||^,  which,  for  our  choice  of  uq  is  the 
same  as  Ua  = 2 ||fc||oo-  This  is  closely  related  to  the  value  used  by  McShane  et  al.  [McSMS88| 
of  Ua  = 1 + 2 ||rj|^.  Mehrotra  [Meh89],  on  the  other  hand,  sets  Ua  = 2 | mini<t<„{ri}  |.  Both 
values  guarantee  initial  feasibility.  To  obtain  the  results  reported  in  §4,  we  set  Ug  using  both  the 
range  of  the  residuals,  and  the  geometric  mean  of  the  maximum  and  minimum  absolute  values 
of  the  residuals.  The  results  reported  in  §4  were  obtained  by  conditionally  specifying  Ua  based 
on  the  size  of  the  ratio  | rmin/^max  1-  Specifically,  we  used 


where 


I ^min  I if  10 X | Tmin  | < | Tmax 

I Tmin  I +\2<P2{r)  Otherwise, 


I ^min  I "(■  1 ^max  i)  ^ 

1P2(^)  ~ max{l,  (roiax  ^min)}i 

and  where  Ai  = 1 and  A2  = .3.  This  value  of  ensures  that  the  starting  values  uq  and 
Ua  are  safely  interior  without  excessively  inflating  the  original  polytope.  Note  that  ipi(r)  is 
the  geometric  mean  of  the  minimum  and  maximum  residual  after  adding  the  minimum  amount 
required  to  make  all  residuals  greater  than  or  equal  to  one,  and  <P2{t’)  is  the  range  of  the  residuals. 


Problem  Preprocessing.  The  Netlib  test  problem  set  [Gay85|  has  test  problems  that  may 
contain  empty  rows  or  columns,  and  variables  that  are  explicitly  fixed.  A small  amount  of  pre- 
processing is  performed  on  the  test  problems  to  remove  these  extraneous  items.  The  process  first 
identifies  fixed  variables  and  removes  them  from  the  problem.  This  is  a multiple  pass  process, 
since  removal  of  a variable  in  an  earlier  iteration  may  result  in  identifying  additional  fixed  vari- 
ables. Once  all  identified  fixed  variables  are  removed,  the  procedure  removes  all  corresponding 
empty  rows  from  the  problem.  (See  also  [BreMW75].) 

Data  Structures.  The  A matrix  is  stored  in  sparse  format  using  the  XMP  experimental 
mathematical  programming  data  structures  described  in  [Mar8l].  The  Hessian  matrix  A^ A 
is  stored  using  the  data  structures  from  the  Yale  Sparse  Matrix  Package  SMPAK  [SCA85]. 
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Hessian  GMW  Modification.  The  Hessian  is  encoded  and  factored  using  the  Yale  Sparse 
Matrix  Package  SMPAK  [SCA85|.  The  minimum  degree  ordering  subroutine  is  invoked  at  the 
beginning  of  Phase  1 and  again  at  the  beginning  of  Phzise  2 to  find  a permutation  of  the  rows 
and  columns  that  reduces  fill-in.  The  symbolic  factorization  routine  is  also  invoked  only  once 
at  the  beginning  of  each  Phase.  As  discussed  below,  constraints  are  implicitly  dropped  <is  our 
algorithm  progresses.  We  do  not  repeat  the  symbolic  factorization  as  constraints  axe  dropped 
however. 

At  each  iteration,  the  Cholesky  factorization  U^EU  is  formed  using  SMPAK  routine  SNF. 
We  have  altered  this  code  to  detect  near  singularity  of  the  Hessian  using  the  modified  Cholesky 
algorithm  specified  on  page  111  of  GiU  et  al.  [GilMWSl).  In  this  algorithm,  the  Cholesky  factors 
U and  E are  computed  subject  to  two  requirements.  First,  all  elements  Ea  must  be  strictly 
positive.  Second,  the  elements  [U’^ E^!"^)  must  satisfy  the  uniform  bound  specified  in  [GilMWSl]. 
As  the  algorithm  is  presented  in  [GUMWSl],  when  these  conditions  cire  not  met  for  a row  z,  the 
diagonal  elements  of  the  original  matrix  are  implicitly  increased  until  the  conditions  eire  satisfied. 
This  results  in  an  increase  in  the  value  for  Ea,  while  the  entries  of  the  corresponding  row  i of 
U are  left  unchanged. 

In  our  implementation,  we  use  the  criteria  of  the  modified  Cholesky  algorithm  to  identify  and 
remove  rows  of  U that  are  effectively  dependent.  Thus,  the  algorithm  was  changed  so  that  when 
a nonzero  quantity  would  be  added  to  the  diagonal  entry  Ea,  the  corresponding  off-diagonal 
entries  in  U are  zeroed  as  well.  Note  that  this  procedure  is  slightly  more  expensive  than  that 
of  simply  chcinging  nonpositive  diagonal  entries  Ea  to  some  small  positive  number  during  the 
factorization.  It  has  the  advantage,  however,  of  isolating  dependent  columns  of  DA  and, 
thus,  obtaining  solutions  unaffected  by  the  dependent  columns. 

Constraint  Dropping.  Constraints  that  are  sufficiently  far  from  the  current  point  u,  i.e., 
those  having  residuals  ry(u)  that  satisfy 

rj  (u)  > 10^^  X min{rfc(u),  k = 1, . . . , m},  (3.2) 

aire  explicitly  removed  from  the  computations.  Constraints  j that  satisfy  (3.2)  are  “dropped” 
by  setting  Rj  and  Djj  to  zero  prior  to  computing  A^ A and  A^ R.  This  can  improve  the 
sparsity  in  A^ A and  the  numerical  accuracy  of  the  resulting  search  directions,  and  therefore 
can  lead  to  improved  performance. 

Steplength  Selection.  All  procedures  reported  in  §4  use  a steplength  equal  to  99%  of  the 
maximum  feasible  steplength  in  the  designated  search  direction.  This  steplength  shifts  the 
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current  estimate  u to  within  1%  of  the  distance  to  the  boundary  of  the  polytope. 


Stopping  Criteria.  Three  standard  convergence  tests  are  used  to  terminate  the  iterations. 
These  tests  are  based  on  (a)  the  relative  change  in  the  objective  function,  (b)  the  relative 
difference  between  the  primal  and  dual  objective  values  (see,  e.g.,  [AdlRV86|),  and  (c)  the 

steplength. 

Objective  function  convergence  is  obtained  when 


If!— < 10-8  (3 

I I 

where  Zi  is  the  objective  function  value  at  iteration  i.  The  convergence  criterion  based  on  the 
relative  difference  between  the  primal  and  dual  objective  values  is 


zf  - Zi 


max{|  zf  I,  I Zi  1} 


< 10 


-8 


(3.4) 


where  zf  is  the  dual  objective  function  value  at  the  current  iteration.  This,  of  course,  can 
only  be  tested  when  the  ducil  multiplier  estimate  D^A  [A'^D'^A)  ^ c is  nonnegative  (see,  e.g., 
[MonM87]).  The  third  convergence  criterion  is  based  on  steplength,  where  convergence  is  ob- 
served when  the  steplength  is  less  than  or  equal  to  10“^®.  All  problems  in  our  test  set  met  either 
(3.3)  or  (3.4). 


Computing  Environment.  The  methods  reported  here  are  implemented  in  Fortran  77  and 
executed  in  double  precision  on  a Sun  3/60  using  compiler  version  F77-Sun  1.1  and  options  -O 
and  -f68881. 


4.  Computational  Results 

During  thk  study,  several  variants  of  the  methods  described  in  §2  were  analyzed.  In  this  section, 
computational  results  are  presented  for  the  optimal  3-dimensional  method  and  for  two  variants  of 
the  dual  affine  method.  Previously,  it  has  been  shown  in  [AdlRV86]^  [MonM87],  and  [McSMS88] 
that  the  dual  affine  method  compares  favorably  to  MINOS  5.0  [MurS83],  a weU-known  and 
widely-available  implementation  of  the  simplex  method.  No  direct  comparison  with  MINOS  is 
made,  however,  in  this  paper. 

The  methods  analyzed  in  this  study  were  tested  on  31  of  the  54  publicly  available  linear 
programming  problems  available  through  Netlib  [Gay85].  The  problems  omitted  from  our  study 
are  those  with  explicit  bounds,  which  our  implementations  do  not  currently  handle.  The  size 
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and  optimal  objective  value  for  each  of  these  31  problems  are  shown  in  Table  4.1.  The  number 
of  rows  and  columns  specified  for  the  Original  Problem,  as  well  as  the  optimal  objective  value, 
are  supplied  by  [Gay85].  These  problem  dimensions  exclude  the  slack  column,  the  right-hand 
side  vector,  and  the  cost  row.  The  Reformulated  Problem  dimensions  indicate  the  size  of  the 
dual  problems  that  we  actually  solved  after  preprocessing  the  data  as  described  in  §3. 

All  but  3 of  the  31  problems  analyzed  require  a Phase  1 procedure  to  obtain  an  initial 
feasible  point  when  using  the  starting  values  specified  in  §3.  Another  8 problems  do  not  have 
a fuU  dimensional  interior  and  therefore  only  require  Phase  1 to  find  the  optimal  solution.  The 
remaining  20  problems  require  both  Phase  1 and  Phase  2. 

1 

4.1.  Results 

The  computational  results  provided  here  compare  our  Optimal  S-D  algorithm  with  our  Dual 
Affine  results.  Both  methods  use  the  same  Fortran  “kernel”  routines.  Although  only  the  essential 
computations  are  executed  for  each  method,  we  expect  that  a further  reduction  in  execution 
time  for  both  methods  could  be  obtained  by  customizing  each.  It  may  not  be  meaningful, 
therefore,  to  compare  directly  these  times  with  times  obtained  using  other  codes.  In  addition, 
we  recognize  that  our  timing  results  are  influenced  by  the  order  of  computations,  the  compiler 
and  the  computing  environment.  Even  the  relative  timing  results  reported  here  may  not  stay 
constant  if  run  on  a different  computer.  For  a detailed  discussion  of  the  difficulties  inherent  in 
comparing  numerical  methods  see  [JacBNP89]. 

Table  4.2  summarizes  our  Original  Dual  Affine  results,  computed  as  discussed  in 
[BogDDW89],  and  our  New  Dual  Affne  results  and  our  Optimal  3-D  results,  computed  as 
described  in  §3.  For  each  of  these  methods.  Table  4.2  reports: 

Iter.:  the  number  of  iterations  required  for  Phase  1 and  the  total  number  of  iterations,  with 

the  two  values  listed  as  p/q-, 

Time:  the  CPU  time,  in  seconds,  excluding  data  input,  preprocessing  and  scaling  as  described 
in  §2,  and  formulation  of  the  dual  problem;  and 

Error:  the  relative  error  of  the  solution,  | [z  — z*)lz*  |,  where  z is  the  estimated  optimal 
objective  value  and  z*  is  the  true  optimal  value  listed  in  Table  4.1. 

The  times  for  the  Original  Dual  Affne  results  cire  labeled  by  0,  the  New  Dual  Affne  results 
by  $,  and  the  Optimal  S-D  results  by  for  consistancy  with  Table  4.3.  The  total  number  of 
iterations  and  the  total  times  required  to  run  all  of  the  problems  are  listed  at  the  bottom  of 
Table  4.2. 
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The  iteration  counts  for  our  Original  Dual  Affine  results  are  nearly  identical  to  the  dual 
affine  results  reported  in  [AdlRV86|,  (MonM87]  and  [McSMS88].  We  have  provided  the  times 
for  this  work  to  allow  our  new  results  to  be  compared  to  our  earlier  work  and  that  reported  by 

others. 

Table  4.3  provides  a comparison  of  our  Optimal  S-D  approach  with  our  Dual  Affine  results. 
For  convenience,  the  CPU  times  for  the  each  of  the  methods  axe  repeated  from  Table  4.2.  In 
addition,  Table  4.3  reports  the  change  in  time  between  the  Optimal  S-D  approach  and  the  Dual 
Affine  approaches,  listed  both  in  seconds  and  as  a relative  percentage  of  the  Dual  Affine  times, 
and  the  ratio  of  the  Dual  Affine  times  to  the  Optimal  S-D  times. 

4.2.  Observations 

Convergence.  Table  4.2  shows  that  our  Optimal  S-D  results  agree  weU  with  the  accepted 
optimal  values  provided  in  [Gay85].  Each  problem  in  the  test  set  is  solved  “correctly”,  converging 
to  the  true  value  with  at  least  6 digits,  and  in  all  but  3 cases  with  8 or  more  digits,  of  agreement. 

Iteration  Counts.  Table  4.2  shows  that  our  Optimal  S-D  approach  results  in  a significant 
reduction  in  the  number  of  iterations  when  compared  to  either  of  our  Dual  Affine  results.  The 
number  of  iterations  required  by  the  Optimal  S-D  method  is  less  than  that  required  by  the  Dual 
Affine  methods  for  all  problems.  The  relative  reductions  range  from  17.6%  to  52.9%  for  the 
individual  problems,  and  over  the  set  of  31  problems  the  total  reduction  in  iterations  is  33.4%. 
Since  our  dual  affine  results  compare  favorably  with  those  reported  in  [AdlRV86],  [MonM87], 
[McSMS88|  and  [LusMS89],  we  believe  that  our  Optimal  S-D  results  are  also  competative  with 
most  current  interior  point  work. 

Timings.  The  Optimal  S-D  approach  results  in  an  overall  reduction  in  CPU  time  of  18.9% 
when  compared  with  our  New  Dual  Affine  results,  and  34.9%  when  compared  with  our  Original 
Dual  Affine  results.  When  compared  with  our  New  Dual  Affine  results,  22  of  the  31  problems 
show  a reduction  in  CPU  time  using  the  Optimal  S-D  approach.  Individual  times  decrease  by 
as  much  as  39.2%,  and  15  of  the  problems  decrease  in  time  by  10%  or  more. 

Of  the  problems  experiencing  an  increase  in  time,  the  Optimal  S-D  approach  causes  a relative 
increaise  of  more  than  10%  for  problems  ScagrZ,  Sharelb,  Scsd6  and  ScsdS.  The  largest  absolute 
increase  is  11.7  seconds  (problem  ScsdS),  which  is  16.5%  of  its  individual  execution  time  (82.8 
seconds)  and  is  only  0.2%  of  the  total  execution  time  for  aU  problems  (6562  seconds). 

It  ^ easily  shown  that  the  solution  of  the  subproblem  does  not  increase  the  order  of  work  per 
iteration.  Our  timings  show  that  the  additional  time  required  to  set  up  and  solve  the  subproblem 
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is  only  13.7%  (896  seconds)  of  the  total  CPU  time  (6562  seconds)  for  the  Optimal  S-D  method. 
Of  this,  640  seconds  are  spent  in  the  simplex  solver,  and  the  remaining  time  is  primarily  that 
required  to  solve  for  32  = D^A)  ^ Aj  or  ^ A^i2  and  33  = Sh- 

Data  input  and  processing  requires  a total  of  about  498  seconds  for  the  31  problems.  Of  this, 
123  seconds  are  required  for  the  scaling  and  preprocessing  described  in  §3.  Since  the  data  input 
and  processing  time  is  constant  for  each  of  our  methods,  it  is  relatively  larger  for  our  Optimal 
S-D  method  than  for  our  New  Dual  Affine  method.  Our  results  would  change  slightly  in  favor 
of  our  New  Dual  Affine  results,  therefore,  if  the  data  input  and  processing  times  were  included. 
The  change  is  not  significant,  however,  and  does  not  alter  our  conclusions. 

\ 

4.3.  Conclusions 

This  study  demonstrates  the  computational  advantages  of  using  optimal  third-order  methods  as 
more  sophisticated  adaptations  to  the  traditional  method  of  centers.  In  particular,  our  Optimal 
S-D  results  are  a significant  improvement  over  our  New  Dual  Affine  results.  The  total  number 
of  iterations  is  reduced  by  33.4%,  and  the  toted  CPU  time  by  18.9%.  The  iteration  counts  for 
our  Optimal  S-D  results  are  also  very  competitive  with  the  dual  affine  and  primal-dual  results 
reported  in  [AdlRV86],  [MonM87],  [McSMS88],  and  (LusMS89j. 
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Table  4.1:  Test  Problem  Characteristics 


Problem 

Name 

Original 

Rows  Col. 

Rows 

— B 

Col. 

.eformulated  Prob 

Non. 

Phase  1 

lem  — 

seros 

Phase  2 

Optimal  Value 

Const. 

Hess. 

Const. 

Hess. 

Afiro 

27 

32 

51 

27 

153 

118 

102 

90 

-.46475314286e3 

ADLittle 

56 

97 

137 

55 

554 

433 

417 

377 

.22549496316e6 

Scagr7 

129 

140 

184 

128 

641 

753 

457 

606 

-.23313892547e7 

Sc205 

205 

203 

315 

203 

978 

858 

663 

654 

-.52202061211e2 

Share2b 

96 

79 

162 

96 

939 

968 

777 

871 

-.41573224074e3 

Sharelb 

117 

225 

248 

112 

1396 

1080 

1148 

967 

-.76589318579e5 

Scorpion 

388 

358 

453 

375 

1913 

2367 

1460 

1991 

.18781248227e4 

Scagr25 

471 

500 

670 

470 

2387 

2841 

1717 

2370 

-.14753433060e8 

ScTapl 

300 

480 

660 

300 

2532 

1987 

1872 

1686 

.14122500000e4 

t BrandY 

220 

249 

247 

137 

2227 

2364 

— 

— 

.15185098965e4 

J Scsdl 

77 

760 

760 

77 

— 

— 

2388 

1133 

.86666666743el 

brael 

174 

142 

316 

174 

2759 

11402 

2443 

11227 

-.89664482186e6 

BandM 

305 

472 

425 

258 

2459 

3370 

2034 

3111 

-.15862801845e3 

t Scfxml 

330 

457 

592 

322 

3299 

3526 

— 

— 

.18416759028e5 

tE226 

223 

282 

469 

220 

3206 

3012 

— 

— 

-.18751929066e2 

t Scrs8 

490 

1169 

1250 

465 

4428 

2450 

— 

— 

.90429695380e3 

t Beaconfd 

173 

262 

262 

140 

3357 

2504 

— 

— 

.33592485807e5 

J Scsd6 

147 

1350 

1350 

147 

— 

— 

4316 

2099 

.50500000078e2 

Ship04s 

402 

1458 

1414 

268 

5538 

3265 

4124 

2996 

.17987147004e7 

t Scfxm2 

660 

914 

1184 

644 

6603 

7071 

— 

— 

.36660261564e5 

Ship04l 

402 

2118 

2162 

356 

8530 

4933 

6368 

4576 

.17933245379e7 

ShipOSs 

778 

2387 

2171 

416 

8477 

4969 

6306 

4552 

.19200982105e7 

ScTap2 

1090 

1880 

2500 

1090 

9834 

7686 

7334 

6595 

.17248071428e4 

t ScfxmS 

990 

1371 

1776 

966 

9907 

10616 

— 

— 

.54901254549e5 

Ship  12s 

1151 

2763 

2293 

466 

8849 

5126 

6556 

4659 

.14892361344e7 

J ScsdS 

397 

2750 

2570 

397 

— 

— 

8584 

4280 

.90499999992e3 

ScTapS 

1480 

2480 

3340 

1480 

13074 

10347 

9734 

8866 

.14240000000e4 

CzProb 

929 

3523 

3141 

737 

12595 

7715 

9454 

6977 

.21851966989e7 

t 25FV47 

821 

1571 

1849 

793 

12415 

12509 

— 

— 

.55018458883e4 

ShipOSl 

778 

4283 

4339 

688 

17149 

9841 

12810 

9152 

.19090552113e7 

Ship  121 

1151 

5427 

5329 

838 

20993 

11942 

15664 

11103 

.14701879193e7 

t indicates  Phcise  1 problem 
t indicates  Phase  2 problem 


Table  4.2:  Results 


Problem 

Original  Dual  Affine 

Iter.  Time  Error 

(sec) 

New  Dual  Affine 

Iter.  Time  Error 

(sec) 

Optimal  3-D 

Iter.  Time  Error 

(sec) 

Afiro 

1/ 

20 

1.3 

5.E-09 

2/ 

19 

1.3 

9.E-09 

1/ 

10 

1.3 

4.E-09 

ADLittle 

1/ 

22 

5.3 

4.E-09 

2/ 

25 

6.1 

9.E-10 

1/ 

15 

6.0 

3.E-11 

Scagr? 

3/ 

24 

8.7 

2.E-07 

3/ 

25 

8.9 

2.E-07 

2/ 

18 

10.6 

2.E-07 

Sc205 

4/ 

27 

15.1 

3.E-09 

6/ 

28 

16.9 

2.E-09 

3/ 

17 

15.7 

l.E-10 

Share2b 

4/ 

30 

13.6 

3.E-09 

3/ 

27 

12.4 

3.E-09 

2/ 

14 

9.9 

5.E-09 

Sharelb 

7/ 

37 

27.0 

9.E-09 

4/ 

34 

23.3 

4.E-09 

2/ 

28 

28.4 

5.E-11 

Scorpion 

5/ 

23 

37.7 

4.E-09 

2/ 

26 

33.2 

5.E-09 

2/ 

17 

33.5 

2.E-12 

Scagr25 

3/ 

28 

50.5 

l.E-09 

3/ 

33 

55.7 

3.E-10 

2/ 

21 

55.8 

l.E-11 

ScTapl 

6/ 

34 

55.5 

8.E-09 

2/ 

35 

49.5 

9.E-09 

1/ 

21 

45.9 

9.E-10 

t BrandY 

37/ 

37 

108.0 

l.E-05 

32/ 

32 

74.6 

6.E-09 

21/ 

21 

58.1 

l.E-10 

t Scsdl 

0/ 

17 

19.5 

5.E-09 

0/ 

17 

19.2 

5.E-09 

0/ 

8 

14.5 

2.E-09 

Israel 

9/ 

40 

479.9 

8.E-09 

4/ 

31 

375.5 

l.E-08 

3/ 

24 

307.0 

l.E-10 

BandM 

7/ 

30 

87.5 

8.E-09 

7/ 

30 

74.2 

9.E-09 

4/ 

21 

65.3 

2.E-10 

t Scfxml 

33/ 

33 

138.6 

l.E-09 

38/ 

38 

154.9 

l.E-09 

23/ 

23 

113.3 

2.E-10 

tE226 

37/ 

37 

130.3 

3.E-08 

32/ 

32 

109.7 

9.E-09 

24/ 

24 

101.0 

5.E-07 

t Scrs8 

52/ 

52 

450.0 

3.E-08 

34/ 

34 

261.7 

6.E-09 

25/ 

25 

224.7 

2.E-10 

t BeaconFD 

25/ 

25 

94.8 

6.E-09 

25/ 

25 

73.9 

8.E-09 

17/ 

17 

59.3 

2.E-09 

t Scsd6 

0/ 

19 

38.5 

5.E-09 

0/ 

19 

38.5 

5.E-09 

0/ 

13 

43.0 

l.E-09 

Ship04s 

5/ 

29 

129.0 

3.E-09 

2/ 

24 

83.8 

7.E-10 

1/ 

15 

75.1 

6.E-09 

t Scfxm2 

37/ 

37 

449.8 

7.E-09 

39/ 

39 

463.9 

6.E-09 

29/ 

29 

387.9 

8.E-09 

Ship04l 

4/ 

29 

201.8 

2.E-09 

2/ 

22 

153.5 

2.E-08 

1/ 

17 

152.3 

3.E-10 

ShipOSs 

5/ 

31 

271.0 

2.E-09 

1/ 

24 

122.5 

8.E-10 

1/ 

16 

132.3 

4.E-09 

ScTap2 

6/ 

32 

498.4 

2.E-09 

2/ 

27 

362.1 

5.E-09 

1/ 

16 

285.0 

3.E-09 

t ScfxmS 

38/ 

38 

906.2 

l.E-09 

40/ 

40 

896.3 

8.E-09 

31/ 

31 

783.5 

3.E-10 

Shipl2s 

5/ 

30 

363.3 

2.E-08 

1/ 

24 

114.0 

3.E-09 

1/ 

16 

122.3 

2.E-09 

t ScsdS 

0/ 

19 

73.6 

8.E-09 

0/ 

19 

71.1 

8.E-09 

0/ 

13 

82.8 

2.E-10 

ScTapS 

6/ 

34 

766.6 

5.E-09 

2/ 

29 

526.6 

6.E-09 

1/ 

18 

445.9 

4.E-09 

CzProb 

3/ 

44 

970.7 

l.E-09 

2/ 

53 

933.1 

9.E-10 

1/ 

41 

852.0 

3.E-11 

t 25FV47 

55/ 

55 

2443.2 

4.E-08 

51/ 

51 

2095.8 

5.E-08 

28/ 

28 

1275.2 

4.E-06 

ShipOSl 

4/ 

28 

494.8 

l.E-09 

1/ 

27 

375.0 

6.E-10 

1/ 

17 

353.1 

3.E-10 

Shipl2l 

5/ 

28 

755.3 

2.E-08 

2/ 

29 

506.5 

6.E-10 

1/ 

17 

421.4 

l.E-08 

AU 

407/969 

10085.5 

344/918 

8093.7 

230/611 

6562.1 

t indicates  Phase  1 problem 
t indicates  Phcise  2 problem 


Table  4.3:  Comparison 


Optimal  3-D 

New  Dual  Affine 

Original  Dual  Affine 

$ 

^ $ 

4> 

0 

^-0 

0 

0/'F 

Problem 

(sec) 

(sec) 

(sec) 

(%) 

(sec) 

(sec) 

{%) 

Afixo 

1.3 

1.3 

0.0 

0.0 

1.00 

1.3 

0.0 

0.0 

1.00 

ADLittle 

6.0 

6.1 

-0.1 

-1.6 

1.02 

5.3 

0.7 

13.2 

0.88 

Scagr? 

10.6 

8.9 

1.7 

19.1 

0.84 

8.7 

1.9 

21.8 

0.82 

Sc205 

15.7 

16.9 

-1.2 

-7.1 

1.08 

15.1 

0.6 

4.0 

0.96 

Share2b 

9.9 

12.4 

-2.5 

-20.2 

1.25 

13.6 

-3.7 

-27.2 

1.37 

Sharelb 

28.4 

23.3 

5.1 

21.9 

0.82 

27.0 

1.4 

5.2 

0.95 

Scorpion 

33.5 

33.2 

0.3 

0.9 

0.99 

37.7 

-4.2 

-11.1 

1.13 

Scagr25 

55.8 

55.7 

0.1 

0.2 

1.00 

50.5 

5.3 

10.5 

0.91 

ScTapl 

45.9 

49.5 

-3.6 

-7.3 

1.08 

55.5 

-9.6 

-17.3 

1.21 

t BrandY 

58.1 

74.6 

-16.5 

-22.1 

1.28 

108.0 

-49.9 

-46.2 

1.86 

1 Scsdl 

14.5 

19.2 

-4.7 

-24.5 

1.32 

19.5 

-5.0 

-25.6 

1.34 

Israel 

307.0 

375.5 

-68.5 

-18.2 

1.22 

479.9 

-172.9 

-36.0 

1.56 

BandM 

65.3 

74.2 

-8.9 

-12.0 

1.14 

87.5 

-22.2 

-25.4 

1.34 

t Scfxml 

113.3 

154.9 

-41.6 

-26.9 

1.37 

138.6 

-25.3 

-18.3 

1.22 

tE226 

101.0 

109.7 

-8.7 

-7.9 

1.09 

130.3 

-29.3 

-22.5 

1.29 

t ScrsS 

224.7 

261.7 

-37.0 

-14.1 

1.16 

450.0 

-225.3 

-50.1 

2.00 

t BeaconFD 

59.3 

73.9 

-14.6 

-19.8 

1.25 

94.8 

-35.5 

-37.4 

1.60 

t Scsd6 

43.0 

38.5 

4.5 

11.7 

0.90 

38.5 

4.5 

11.7 

0.90 

Ship04s 

75.1 

83.8 

-8.7 

-10.4 

1.12 

129.0 

-53.9 

-41.8 

1.72 

t Scfxm2 

387.9 

463.9 

-76.0 

-16.4 

1.20 

449.8 

-61.9 

-13.8 

1.16 

Ship04l 

152.3 

153.5 

-1.2 

-0.8 

1.01 

201.8 

-49.5 

-24.5 

1.33 

ShipOSs 

132.3 

122.5 

9.8 

8.0 

0.93 

271.0 

-138.7 

-51.2 

2.05 

ScTap2 

285.0 

362.1 

-77.1 

-21.3 

1.27 

498.4 

-213.4 

-42.8 

1.75 

t ScfxmS 

783.5 

896.3 

-112.8 

-12.6 

1.14 

’ 906.2 

-122.7 

-13.5 

1.16 

Ship  12s 

122.3 

114.0 

8.3 

7.3 

0.93 

363.3 

-241.0 

-66.3 

2.97 

1 ScsdS 

82.8 

71.1 

11.7 

16.5 

0.86 

73.6 

9.2 

12.5 

0.89 

ScTapS 

445.9 

526.6 

-80.7 

-15.3 

1.18 

766.6 

-320.7 

-41.8 

1.72 

CzProb 

852.0 

933.1 

-81.1 

-8.7 

1.10 

970.7 

-118.7 

-12.2 

1.14 

t 25FV47 

1275.2 

2095.8 

-820.6 

-39.2 

1.64 

2443.2 

-1168.0 

-47.8 

1.92 

ShipOSl 

353.1 

375.0 

-21.9 

-5.8 

1.06 

494.8 

-141.7 

-28.6 

1.40 

Ship  121 

421.4 

506.5 

-85.1 

-16.8 

1.20 

755.3 

-333.9 

-44.2 

1.79 

AU 

6562.1 

8093.7 

-1531.6 

-18.9 

1.23 

10085.5 

-3523.4 

-34.9 

1.54 

t indicates  Phase  1 problem 
J indicates  Phase  2 problem 
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